## Simulation Flags
RunMCMC <- FALSE
MCMCBurn <- 1000
MCMCSample <- 4000
MCMCThin <- 2
MCMCChains <- 2

#library(eliot)
library(Rcpp)
library(inline)
library(rjags)
library(car)
library(KernSmooth)
load.module('dic')
load(file='Output/VoteDetailsForFittedSubset.Rdata')
load(file='Output/CARIRTSavedOutput.Rdata')

load('Data/SCDB_2012_01_justiceCentered_Citation.RData')
SCDB <- SCDB_2012_01_justiceCentered_Citation


pdf("Output/ConvergenceChecks.pdf",12,8)
par(mfrow=c(2,2))

plot(carirt.out$lambda.chain[1,],type="l",main="Lambda",xlab="Iteration",ylab="")

plot(0,0,type="n",xlim=c(0,length(carirt.out$lambda.chain)+1),ylim=c(-0.75,0.75),xlab="Iteration",ylab="",main="Average Justice Preferences")
for (i in 1:32) lines(1:length(carirt.out$lambda.chain),carirt.out$average.pref.chain[i,],type="l",col=rgb(0,0,0,0.1))

# plot(0,0,type="n",xlim=c(0,length(carirt.out$lambda.chain)+1),ylim=c(-0.75,0.75),xlab="Iteration",ylab="",main="Scalia's Preferences in Kyllo v US")
# lines(1:length(carirt.out$lambda.chain),carirt.out$latent.pref.chain[6858,27,],type="l",col=rgb(0,0,0,0.1))

plot(0,0,type="n",xlim=c(0,length(carirt.out$lambda.chain)+1),ylim=c(-1,1),xlab="Iteration",ylab="",main="Polarity of Kyllo v US")
lines(1:length(carirt.out$lambda.chain),carirt.out$beta.chain[6858,],type="l",col=rgb(0,0,0,0.1))

plot(0,0,type="n",xlim=c(0,length(carirt.out$lambda.chain)+1),ylim=c(-1,1),xlab="Iteration",ylab="",main="Cutpoint of Kyllo v US")
lines(1:length(carirt.out$lambda.chain),carirt.out$alpha.chain[6858,],type="l",col=rgb(0,0,0,0.1))

dev.off()


#############################################
## Fill out additional variables from SCDB ##
#############################################

n.cases <- dim(carirt.out$latent.pref.chain)[1]
n.voters <- dim(carirt.out$latent.pref.chain)[2]
n.sims <- dim(carirt.out$latent.pref.chain)[3]
voter.names <- dimnames(carirt.out$latent.pref.chain)[[2]]

colnames(Y) <- voter.names

voteDetails$Chief <- rep("", n.cases)
voteDetails$IssueArea <- rep(NA, n.cases)
voteDetails$Term <- rep(NA, n.cases)
voteDetails$Author <- rep(NA, n.cases)
voteDetails$Assigner <- rep(NA, n.cases)
voteDetails$MajorityVotes <- rep(NA, n.cases)

JoinMatrix <- matrix(NA,nrow= n.cases, ncol= n.voters)
colnames(JoinMatrix) <- voter.names
rownames(JoinMatrix) <- voteDetails$caseId
SCDB$join <- recode(SCDB$vote,"1=1;2=NA;3=1;4=0;5=1;6:9=NA")

for(i in 1:n.cases){
	find.scdb.rows <- which(voteDetails$caseId[i] == SCDB$caseId)
	voteDetails$IssueArea[i] <- SCDB$issueArea[find.scdb.rows[1]]	
	voteDetails$Term[i] <- SCDB$term[find.scdb.rows[1]]		
	voteDetails$Chief[i] <- SCDB$chief[find.scdb.rows[1]]	
	voteDetails$Author[i] <- SCDB$majOpinWriter[find.scdb.rows[1]]
	voteDetails$Assigner[i] <- SCDB$majOpinAssigner[find.scdb.rows[1]]
	voteDetails$MajorityVotes[i] <- SCDB$majVotes[find.scdb.rows[1]]	
	for (ii in find.scdb.rows){
		tempJustice <- which(SCDB[ii,]$justiceName == colnames(JoinMatrix) )
		JoinMatrix[i,tempJustice] <- SCDB$join[ii]
	}		
}


IssueAreaNames <- c("Criminal Procedure","Civil Rights","First Amendment","Due Process","Privacy","Attorneys","Unions","Economic Activity","Judicial Power","Federalism","Interstate Relations","Federal Taxation","Miscellaneous")
voteDetails$IssueArea <- factor(voteDetails$IssueArea,labels= IssueAreaNames)

voteDetails$Chief <- replace(voteDetails$Chief, voteDetails$Chief == "Vinson","FMVinson")
voteDetails$Chief <- replace(voteDetails$Chief, voteDetails$Chief == "Warren","EWarren")
voteDetails$Chief <- replace(voteDetails$Chief, voteDetails$Chief == "Burger","WEBurger")
voteDetails$Chief <- replace(voteDetails$Chief, voteDetails$Chief == "Rehnquist","WHRehnquist")

voteDetails$Chief <- factor(voteDetails$Chief,levels=colnames(Y))

justiceNameIndex <- rep("",120)
counter <- 1
for (i in unique(SCDB$justice)){
		justiceNameIndex[i] <- unique(SCDB$justiceName)[counter]
		counter <- counter + 1
	}

voteDetails$Author <- justiceNameIndex[voteDetails$Author]
voteDetails$Assigner  <- justiceNameIndex[voteDetails$Assigner]
voteDetails$Author <- factor(voteDetails$Author,colnames(Y))
voteDetails$Assigner  <- factor(voteDetails$Assigner,colnames(Y))

voteDetails$NaturalCourt <- factor(voteDetails$NaturalCourt)
voteDetails$ChiefID <- as.factor(trunc(as.numeric(as.character(voteDetails$NaturalCourt))/100))

ChiefOfNaturalCourt <- rep(NA,length(levels(voteDetails$NaturalCourt )))
for (i in 1:length(ChiefOfNaturalCourt)) {
	ChiefOfNaturalCourt[i] <- as.numeric(voteDetails$ChiefID)[which(as.numeric(voteDetails$NaturalCourt) == i)[1]]
	}
	
latent.preferences <- apply(carirt.out$latent.pref.chain,c(1,2),mean)
cutpoints.est <- apply(carirt.out$alpha.chain,c(1),mean)
polarities.est <- apply(carirt.out$beta.chain,c(1),mean)	

onCourt <- replace(Y,Y == 0,1)
NotOnCourt <- is.na(onCourt)
inMajority <- replace(Y,Y == 0,NA)
NotInMajority <- is.na(inMajority)*1

MajorityVotes <- voteDetails$MajorityVotes
Author <- as.numeric(voteDetails$Author)

ChiefInMajority <- rep(NA,dim(voteDetails)[1])
for(i in 1:dim(voteDetails)[1]){
	ChiefInMajority[i] <- Y[i,as.numeric(voteDetails$Chief)[i]]	
}
ChiefInMajority <- replace(ChiefInMajority,is.na(ChiefInMajority),0)

## clean up justice names 
voter.names <- gsub("[A-Z]+([A-Z])","\\1",voter.names)
voter.names[voter.names == "Connor"] <- "O'Connor"
voter.names[voter.names == "Harlan2"] <- "Harlan"



## Plot average justice positions over time...

pdf('Output/TimeTrendPlot.pdf',8,4.8)

highlight <- is.element(1:34,c(1,17,22,29))

plot(0,0,type="n",ylim=c(-1,1),xlim=c(1946,2010),xlab="Year",ylab="Average Latent Preferences",axes=FALSE,main="Average Case-Specific Preferences")
axis(1,at=seq(1950,2010,10))
axis(2)
for (i in 1:32){
	keepvotes <- which(!is.na(onCourt[,i]))
	tempprefs <- latent.preferences[keepvotes,i]
	#ks.out <- ksmooth(voteDetails$Term[keepvotes],tempprefs,kernel="normal",bandwidth=4)
	#lines(ks.out,col=rgb(0,0,0,0.1+0.4*highlight[i]),lwd=2)
	#text(rev(ks.out$x)[1],rev(ks.out$y)[1], voter.names[i],pos=4,cex=0.5,col=rgb(0,0,0,0.5+0.5*highlight[i]))
	locpoly.out <- locpoly(voteDetails$Term[keepvotes],tempprefs,degree=1,kernel="normal",bandwidth=2)
	lines(locpoly.out,col=rgb(0,0,0,0.1+0.4*highlight[i]),lwd=2)
	text(rev(locpoly.out$x)[1],rev(locpoly.out$y)[1], voter.names[i],pos=4,cex=0.5,col=rgb(0,0,0,0.5+0.5*highlight[i]))
	}

dev.off()

## Plot MQ score equivalents...

pdf('Output/TimeTrendPlotMQ.pdf',8,4.8)

MQ <- read.csv("Data/justices.csv")
MQ <- subset(MQ,is.element(MQ$term,1946:2004))

highlightset <- c("White","Black","Blackmun","Souter")

plot(0,0,type="n",ylim=c(-7,7),xlim=c(1946,2010),xlab="Year",ylab="Martin-Quinn Scores",axes=FALSE,main="Martin-Quinn Scores")
axis(1,at=seq(1950,2010,10))
axis(2)
for (i in unique(MQ$justiceName)){
	tmpterm <- MQ$term[i == MQ$justiceName]
	tmpscore <- MQ$post_mn[i == MQ$justiceName]
	locpoly.out <- locpoly(tmpterm, tmpscore,degree=1,kernel="normal",bandwidth=2)
	lines(locpoly.out,col=rgb(0,0,0,0.1+0.4*(is.element(i,highlightset))),lwd=2)
	text(rev(locpoly.out$x)[1],rev(locpoly.out$y)[1], i,pos=4,cex=0.5,col=rgb(0,0,0,0.5+0.5*(is.element(i,highlightset))))
	
	#lines(tmpterm,tmpscore,col=rgb(0,0,0,0.1+0.4*(is.element(i,highlightset))),lwd=2)
	#text(rev(tmpterm)[1],rev(tmpscore)[1],i,pos=4,cex=0.5,col=rgb(0,0,0,0.5+0.5*is.element(i,highlightset)))
	}

dev.off()





######################################################################
## Calculate Court Median, Majority Median, and Chief's Preferences ##
######################################################################

### Calculate posterior means of distances

MedianDistance.chain <- array(NA,dim=c(n.cases,n.voters,n.sims))
MajMedianDistance.chain <- array(NA,dim=c(n.cases,n.voters,n.sims))
ChiefDistance.chain <- array(NA,dim=c(n.cases,n.voters,n.sims))
CutDistance.chain <- array(NA,dim=c(n.cases,n.voters,n.sims))
AuthorDistance.chain <- array(NA,dim=c(n.cases,n.voters,n.sims))

for (sim in 1:n.sims){
	
	print(paste("Constructing quantities of interest for iteration",sim,"of",n.sims,"..."))
	
	latent.pref.sim <- carirt.out$latent.pref.chain[,,sim]

	# Median Distance

	court.median.sim <- apply(latent.pref.sim * onCourt,1,median,na.rm=TRUE)
	MedianDistance.chain[,,sim] <- abs(latent.pref.sim - matrix(court.median.sim,nrow=dim(Y)[1],ncol=dim(Y)[2]))

	# Majority Median Distance

	majority.median.sim <- apply(latent.pref.sim* inMajority,1,median,na.rm=TRUE)
	MajMedianDistance.chain[,,sim] <- abs(latent.pref.sim - matrix(majority.median.sim,nrow=dim(Y)[1],ncol=dim(Y)[2]))

	# Chief Distance

	extractmat <- cbind(x=1:n.cases,y=as.numeric(voteDetails$Chief))
	chief.preferences.sim <- latent.pref.sim[extractmat]
	ChiefDistance.chain[,,sim] <- abs(latent.pref.sim - matrix(chief.preferences.sim,nrow=dim(Y)[1],ncol=dim(Y)[2]))

	# Cutline Distance

	cutpoint.sim <- carirt.out$alpha.chain[,sim]
	CutDistance.chain[,,sim] <- abs(latent.pref.sim - matrix(cutpoint.sim,nrow=dim(Y)[1],ncol=dim(Y)[2]))

	# Author Distance
	
	extractmat <- cbind(x=1:n.cases,y=Author)
	author.preferences.sim <- latent.pref.sim[extractmat]
	AuthorDistance.chain[,,sim] <- abs(latent.pref.sim - matrix(author.preferences.sim,nrow=dim(Y)[1],ncol=dim(Y)[2]))

}

MedianDistance <- apply(MedianDistance.chain,c(1,2),mean)
MajMedianDistance <- apply(MajMedianDistance.chain,c(1,2),mean)
ChiefDistance <- apply(ChiefDistance.chain,c(1,2),mean)
CutDistance <- apply(CutDistance.chain,c(1,2),mean)
AuthorDistance <- apply(AuthorDistance.chain,c(1,2),mean)

# add small value to avoid infinite density in measurement error model for chief and author distance measures
MedianDistanceSE <- apply(MedianDistance.chain,c(1,2),sd) + 0.001
MajMedianDistanceSE <- apply(MajMedianDistance.chain,c(1,2),sd) + 0.001
ChiefDistanceSE <- apply(ChiefDistance.chain,c(1,2),sd) + 0.001
CutDistanceSE <- apply(CutDistance.chain,c(1,2),sd) + 0.001
AuthorDistanceSE <- apply(AuthorDistance.chain,c(1,2),sd) + 0.001



court.median.est <- apply(latent.preferences* onCourt,1,median,na.rm=TRUE)




for (j in 1:n.voters){

labelside <- sign(mean(latent.preferences[which(!is.na(onCourt[,j])),j]))

pdf(paste0("Output/Justice", voter.names[j],"ByIssueArea.pdf"),6,5)

plot(0,0,type="n",ylim=c(0,14),xlim=c(-2,2),xlab="Latent Preferences",ylab="Issue Area",main=paste("Justice",colnames(onCourt)[j],"by Issue Area"),axes=FALSE)
axis(1)
for (i in 1:13){
	keepvotes <- which(!is.na(onCourt[,j]) & as.numeric(voteDetails$IssueArea) == i)
	tempprefs <- latent.preferences[keepvotes,j]
	points(tempprefs,rep(i,length(tempprefs)),pch="|",col=rgb(1,0,1,0.05),cex=0.75)
	points(mean(tempprefs),i,pch=16,col=rgb(0,0,0,0.75),cex=0.75)	
	text(-2* labelside,i, levels(voteDetails$IssueArea)[i] ,pos=3+labelside)
	}

dev.off()

}

fadecol <- function(col,alpha){
	components <- col2rgb(col)
	return(rgb(components[1,]/255,components[2,]/255,components[3,]/255,alpha=alpha))
	}	

pdf('Output/IssueMedianTrendPlot.pdf',5,4)

plot(0,0,type="n",ylim=c(-0.5,0.5),xlim=c(0,length(MajorityVotes)+1000),xlab="Case, Sorted by Decision Order",ylab="Average Latent Preferences")
for (i in 1:12){
	keepvotes <- which(as.numeric(voteDetails$IssueArea) == i)
	tempprefs <- court.median.est[keepvotes]
	ks.out <- ksmooth(keepvotes,tempprefs,kernel="normal",bandwidth=600)
	lines(ks.out,lwd=2,col=i)
	text(rev(ks.out$x)[1]-100,rev(ks.out$y)[1],IssueAreaNames[i],pos=4,cex=0.5,col=i)
	}

dev.off()










############################
## Create Pivotality Plot ##
############################	

pivotal.probs <- matrix(NA,dim(Y)[1],dim(Y)[2])

for (j in 1:dim(Y)[1]){
	if(sum(!is.na(Y[j,])) == 9) pivotal.probs[j,] <- rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 5)
	if(sum(!is.na(Y[j,])) == 8) pivotal.probs[j,] <- (rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 5) + rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 4))/2
	if(sum(!is.na(Y[j,])) == 7) pivotal.probs[j,] <- rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 4)
	if(sum(!is.na(Y[j,])) == 6) pivotal.probs[j,] <- (rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 4) + rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 3))/2	
	if(sum(!is.na(Y[j,])) == 5) pivotal.probs[j,] <- rowMeans(apply(carirt.out$latent.pref.chain[j,,] * matrix(replace(Y[j,],Y[j,] == 0,1),dim(carirt.out$latent.pref.chain)[2],dim(carirt.out$latent.pref.chain)[3]),2,rank) == 3)	
}

pivotal.probs <- replace(pivotal.probs,is.na(Y),NA)

pdf("Output/PivotalProbPlot1.pdf",9,4)
plot(0,0,type="n",xlim=c(1946,2005),ylim=c(0,0.25),xlab="Year",ylab="Fraction Pivotal",main="Fraction of Cases Where Each Justice Was Pivotal")
for (i in 1:dim(Y)[2]){
	iobs <- which(!is.na(Y[,i]))
	saved.smooth <- ksmooth(voteDetails$Term[iobs],pivotal.probs[iobs,i],bandwidth=4,kernel="normal")
	text(saved.smooth$x[1],saved.smooth$y[1],colnames(Y)[i],cex=0.5,col=i)
	text(tail(saved.smooth$x,1),tail(saved.smooth$y,1), voter.names[i],cex=0.5,col=i)
	lines(saved.smooth$x,saved.smooth$y,col=i)
}
dev.off()



pdf("Output/PivotalProbPlot2.pdf",9.5,7)
inflation <- 2
bw <- 3
plot(0,0,type="n",xlim=c(1940,2010),ylim=c(0,dim(Y)[2]+1),xlab="Year",ylab="",main="Fraction of Cases Where Each Justice Was Pivotal",axes=FALSE)
axis(1)
for (i in 1:dim(Y)[2]){
	iobs <- which(!is.na(Y[,i]))
	saved.smooth <- ksmooth(voteDetails$Term[iobs],pivotal.probs[iobs,i],bandwidth= bw,kernel="normal")
	text(saved.smooth$x[1],i, voter.names[i],cex=0.75,pos=2)
	polygon(c(saved.smooth$x,rev(saved.smooth$x)),c(i + inflation*saved.smooth$y,rev(i - inflation *saved.smooth$y)),col=1)
}
dev.off()




pdf("Output/PivotalProbPlot3.pdf",9.5,7)
inflation <- 4
plot(0,0,type="n",xlim=c(1940,2012),ylim=c(0,dim(Y)[2]+1),xlab="Term",ylab="",main="Fraction of Cases Where Each Justice Was Pivotal",axes=FALSE)
axis(1)
for (i in 1:dim(Y)[2]){
	for (yr in 1946:2006){
		iobs <- which((!is.na(Y[,i])) & (voteDetails$Term == yr))
		tempfrac <- mean(pivotal.probs[iobs,i])
		if (!is.na(tempfrac )) rect(yr-0.5,i-0.5,yr+0.5,i+0.5,col=grey((1-inflation*tempfrac)),border=NA)
	}
	text(min(voteDetails$Term[!is.na(Y[,i])]),i, voter.names[i],cex=0.65,pos=2)	
}


for (i in 1:5){
	rect(2006-0.5,i-0.5, 2006 +0.5,i+0.5,col=grey((1-inflation*0.05*(i-1))),border=NA)
	text(2006,i,0.05*(i-1),cex=0.65,pos=2)
}
text(2006,6,"Key",cex=0.65)

dev.off()





pdf("Output/PivotalProbPlot4.pdf",9.5,7)
plot(0,0,type="n",xlim=c(1940,2010),ylim=c(0,1),xlab="Year",ylab="Fraction Pivotal",main="Fraction of Cases Where Each Justice Was Pivotal",axes=FALSE)
axis(1)
fracmatrix <- matrix(NA,dim(Y)[2],2004-1945)
for (i in 1:dim(Y)[2]){
	for (yr in 1946:2004){
		iobs <- which((!is.na(Y[,i])) & (voteDetails$Term == yr))
		fracmatrix[i,yr-1945] <- mean(pivotal.probs[iobs,i])
	}
}	

fracmatrix <- replace(fracmatrix,is.na(fracmatrix),0)	
fracmatrix <- fracmatrix / matrix(colSums(fracmatrix),dim(Y)[2],2004-1945,byrow=TRUE)

ideosort <- order(apply(carirt.out$latent.pref.chain,2,mean))

fracmatrix <- fracmatrix[ideosort,]

cumsummatrix <- apply(fracmatrix,2,cumsum)

for (i in 1:dim(Y)[2]){
	lines(1946:2004, cumsummatrix[i,])
}	
	#if (!is.na(tempfrac )) rect(yr-0.5,i-0.5,yr+0.5,i+0.5,col=grey((1-inflation*tempfrac)),border=NA)
	#text(min(voteDetails$Term[!is.na(Y[,i])]),i, voter.names[i],cex=0.75,pos=2)	

dev.off()










################
## Case Plots ##
################

	
plotvote <- function(selectedvote,title){
	presentvoters <- which(is.na(Y[selectedvote,]) == FALSE)
	plot(0,0,type="n",axes=FALSE,xlim=c(-3.5,3.5),ylim=c(0,10),xlab="Latent Preferences Relative to Cutpoint",ylab="",main=title,cex.main=1)
	axis(1)
	abline(v = 0)	
	for (i in 1:length(presentvoters)){
		temp.mcmc <- (carirt.out$latent.pref.chain[selectedvote, presentvoters[i],] - carirt.out$alpha.chain[selectedvote,])
		temp.quantile <- quantile(temp.mcmc,c(0.025,0.25,0.5,0.75,0.975))
		lines(temp.quantile[c(2,4)],rep(i,2),lwd=2)
		lines(temp.quantile[c(1,5)],rep(i,2),lwd=1)
		points(temp.quantile[3],i,pch=16,cex=0.75)
		if (mean(sign(temp.mcmc)) > 0) text(temp.quantile[5],i,names(presentvoters)[i],pos=4,cex=0.75) 
		if (mean(sign(temp.mcmc)) < 0) text(temp.quantile[1],i,names(presentvoters)[i],pos=2,cex=0.75)  
	}
}	

pdf("Output/CasePlots.pdf",6,9)	
par(mfrow=c(3,2))
plotvote(3088,"Roe v. Wade") # Roe
plotvote(6108,"Planned Parenthood v. Casey") # Casey
plotvote(2173,"Miranda v. Arizona") # Miranda
plotvote(2346,"Katz v. United States") # Katz
plotvote(3198,"Miller v. California") # Miller v. Calif
plotvote(6858,"Kyllo v. United States") # Kyllo
dev.off()

pdf("Output/CasePlots.pdf",7,8)	
par(mfrow=c(2,2))
plotvote(3088,"Roe v. Wade") # Roe
plotvote(6108,"Planned Parenthood v. Casey") # Casey
plotvote(2173,"Miranda v. Arizona") # Miranda
plotvote(6858,"Kyllo v. United States") # Kyllo
dev.off()



pdf("Output/CasePlotRoe.pdf",6,4)	
par(mfrow=c(1,1))
plotvote(3088,"Roe v. Wade") # Roe
dev.off()


pdf("Output/CasePlotCasey.pdf",6,4)	
par(mfrow=c(1,1))
plotvote(6108,"Planned Parenthood v. Casey") # Casey
dev.off()

pdf("Output/CasePlotKatz.pdf",6,4)	
par(mfrow=c(1,1))
plotvote(2346,"Katz v. United States") # Katz
dev.off()

pdf("Output/CasePlotMiranda.pdf",6,4)	
par(mfrow=c(1,1))
plotvote(2173,"Miranda v. Arizona") # Miranda
dev.off()

pdf("Output/CasePlotKyllo.pdf",6,4)	
par(mfrow=c(1,1))
plotvote(6858,"Kyllo v. United States") # Kyllo
dev.off()

pdf("Output/CasePlotMiller.pdf",6,4)	
par(mfrow=c(1,1))
plotvote(3198,"Miller v. California") # Miller
dev.off()










## Free up memory ##
rm(carirt.out,MedianDistance.chain, MajMedianDistance.chain, ChiefDistance.chain, CutDistance.chain, AuthorDistance.chain)
gc()


#####################################
## Generate Descriptive Statistics ##
#####################################

quantile(MedianDistance[!NotInMajority],c(0.25,0.5,0.75))
quantile(MajMedianDistance[!NotInMajority],c(0.25,0.5,0.75))
quantile(ChiefDistance[!NotInMajority],c(0.25,0.5,0.75))
quantile(CutDistance[!NotInMajority],c(0.25,0.5,0.75))


sd(MedianDistance[!NotInMajority])
sd(MajMedianDistance[!NotInMajority])
sd(ChiefDistance[!NotInMajority])
sd(CutDistance[!NotInMajority])




####################################
## Run JAGS Models for Authorship ##
####################################

keepCases <- ifelse(rowSums(is.na(MedianDistance)) + rowSums(is.na(MajMedianDistance)) + rowSums(is.na(ChiefDistance))==0,1,0) & 
	(ChiefInMajority) & 
	(apply(latent.preferences,1,sd) < 2) # removes disjoint cases with preference estimates that are not identified relative to the others

# load data and run simulations
forJags <- list(
	Author = Author[keepCases],
	MedianDistanceEst = MedianDistance[keepCases,],
	MajMedianDistanceEst = MajMedianDistance[keepCases,],
	CutDistanceEst = CutDistance[keepCases,],	
	ChiefDistanceEst = ChiefDistance[keepCases,],
	MedianDistancePrec = MedianDistanceSE[keepCases,]^(-2),
	MajMedianDistancePrec = MajMedianDistanceSE[keepCases,]^(-2),
	CutDistancePrec = CutDistanceSE[keepCases,]^(-2),	
	ChiefDistancePrec = ChiefDistanceSE[keepCases,]^(-2),
	NotInMajority = NotInMajority[keepCases,],
	NaturalCourt = as.numeric(voteDetails$NaturalCourt)[keepCases],
	ChiefOfNaturalCourt = ChiefOfNaturalCourt,
	Ncases = table(keepCases)[2],
	Njustices = dim(MedianDistance)[2],
	Ncourts = length(levels(voteDetails$NaturalCourt)),
	Nchiefs = length(unique(ChiefOfNaturalCourt))
)

ModelNames <- c("Main","1","2","3","4")

if (RunMCMC){
	for (current.model in ModelNames){
	
		current.model.string <- paste("AuthorshipMeasurementError", current.model,".jags",sep="")
	
		model.jags <- jags.model(file=current.model.string
			,data=forJags
			,n.chains=MCMCChains,n.adapt=round(MCMCBurn/4))
		
		update(model.jags, round(3*MCMCBurn/4))
		
		posterior <- jags.samples(model.jags, 
			c("beta","beta2","beta3"
			,"alpha","alpha2","alpha3"
			,"sigmaalpha","sigmaalpha2"
			,"sigmabeta","sigmabeta2"
			,"deviance","pD"
			), 
			n.iter=MCMCSample, thin = MCMCThin, progress.bar="text")
		
		current.posterior.string <- paste("Output/AuthorshipMeasurementError", current.model,".RData",sep="")
			
		save(posterior, file=current.posterior.string)	
		
		DIC <- mean(posterior$deviance) + mean(posterior$pD)	
		print(paste("Model",current.model,"DIC:",DIC))
			
	}	
}


########################
## Analyze Full Model ##
########################	

current.posterior.string <- paste("Output/AuthorshipModelMain.RData",sep="")
load(file=current.posterior.string)	
	
DIC <- mean(posterior$deviance) + mean(posterior$pD)	
print(paste("DIC:",DIC))

alpha3.est <- apply(posterior$alpha3,1,mean)

alpha2.est <- apply(posterior$alpha2,c(1,2),mean)
rownames(alpha2.est) <- colnames(Y)
colnames(alpha2.est) <- as.character(unique(voteDetails$Chief))

alpha.est <- apply(posterior$alpha,c(1,2),mean)
rownames(alpha.est) <- colnames(Y)
colnames(alpha.est) <- as.character(unique(voteDetails$NaturalCourt))

beta3.est <- apply(posterior$beta3,1,mean)
beta3.quantiles <- t(apply(posterior$beta3,1,quantile,c(0.025,0.975)))
print(cbind(beta3.est,beta3.quantiles))

beta2.est <- apply(posterior$beta2,c(1,2),mean)
rownames(beta2.est) <- c("Cutpoint","Court Median","Majority Median","Chief Justice")
colnames(beta2.est) <- as.character(unique(voteDetails$Chief))
#beta2.quantiles <- t(apply(posterior$beta2,c(1,2),quantile,c(0.025,0.975)))
print(beta2.est )

beta.est <- apply(posterior$beta,c(1,2),mean)
rownames(beta.est) <- c("Distance to Cutpoint","Distance to Court Median","Distance to Majority Median","Distance to Chief Justice")
colnames(beta.est) <- as.character(unique(voteDetails$NaturalCourt))
print(beta.est)


#####################
## Plot Full Model ##
#####################	

pdf("Output/AuthorshipCoefficientPlotMultivariate.pdf",10,4)
beta.quantiles <- apply(posterior$beta,c(1,2),quantile,c(0.025,0.975))
n.periods <- dim(beta.est)[2]
par(mfrow=c(1,4),mar=c(4,4,4,4))
for (v in 1:4){
	plot(0,0,xlim=c(-5,5),ylim=c(0,n.periods+4),type="n",axes=FALSE,xlab="Coefficient",main=rownames(beta.est)[v],ylab="Natural Court")
	axis(1)
	abline(v=0)
	for (k in 1:n.periods){
		lines(beta.quantiles[1:2,v,k],rep(k,2),col="grey")
		points(beta.est[v,k],k,pch=16,cex=0.75)
	}
	lines(beta3.quantiles[v,1:2],rep(k+2,2),col="grey",lwd=2)
	points(beta3.est[v],k+2,pch=16,cex=1)


	
	# Chief Justice Divisions #
	text(3,1.5,"Vinson",pos=2,cex=0.75)
	abline(h=2.5)
	text(3,8,"Warren",pos=2,cex=0.75)
	abline(h=13.5)
	text(3,17,"Burger",pos=2,cex=0.75)
	abline(h=20.5)
	text(3,24,"Rehnquist",pos=2,cex=0.75)
	abline(h=27.5)
	text(3,29,"Average",pos=2,cex=0.75)
	
}
	box("outer",col="grey")
dev.off()


############################
## Plot Univariate Models ##
############################	


pdf("Output/AuthorshipCoefficientPlotUnivariate.pdf",10,4)
beta.quantiles <- apply(posterior$beta,c(1,2),quantile,c(0.025,0.975))
n.periods <- dim(beta.est)[2]
par(mfrow=c(1,4),mar=c(4,4,4,4))
for (v in 1:4){
	
	current.posterior.string <- paste("Output/AuthorshipMeasurementError", v,".RData",sep="")
	load(file=current.posterior.string)	
	
	DIC <- mean(posterior$deviance) + mean(posterior$pD)	
	print(paste("DIC:",DIC))
	
	alpha3.est <- apply(posterior$alpha3,1,mean)
	
	alpha2.est <- apply(posterior$alpha2,c(1,2),mean)
	rownames(alpha2.est) <- colnames(Y)
	colnames(alpha2.est) <- as.character(unique(voteDetails$Chief))
	
	alpha.est <- apply(posterior$alpha,c(1,2),mean)
	rownames(alpha.est) <- colnames(Y)
	colnames(alpha.est) <- as.character(unique(voteDetails$NaturalCourt))
	
	beta3.est <- apply(posterior$beta3,1,mean)
	beta3.quantiles <- t(apply(posterior$beta3,1,quantile,c(0.025,0.975)))
	print(cbind(beta3.est,beta3.quantiles))
	
	beta2.est <- apply(posterior$beta2,c(1,2),mean)
	rownames(beta2.est) <- c("Cutpoint","Court Median","Majority Median","Chief Justice")
	colnames(beta2.est) <- as.character(unique(voteDetails$Chief))
	#beta2.quantiles <- t(apply(posterior$beta2,c(1,2),quantile,c(0.025,0.975)))
	print(beta2.est )
	
	beta.est <- apply(posterior$beta,c(1,2),mean)
	beta.quantiles <- apply(posterior$beta,c(1,2),quantile,c(0.025,0.975))
	rownames(beta.est) <- c("Distance to Cutpoint","Distance to Court Median","Distance to Majority Median","Distance to Chief Justice")
	colnames(beta.est) <- as.character(unique(voteDetails$NaturalCourt))
	print(beta.est)

	# Calculateaverage of natural-court specific coefficients
	beta.average.chain <- apply(posterior$beta,c(1,3),mean)
	beta.average.est <- apply(beta.average.chain,1,mean)
	beta.average.quantiles <- apply(beta.average.chain,1,quantile,c(0.025,0.975))
	
	plot(0,0,xlim=c(-5,5),ylim=c(0,n.periods+4),type="n",axes=FALSE,xlab="Coefficient",main=rownames(beta.est)[v],ylab="Natural Court")
	axis(1)
	abline(v=0)
	for (k in 1:n.periods){
		lines(beta.quantiles[1:2,v,k],rep(k,2),col="grey")
		points(beta.est[v,k],k,pch=16,cex=0.75)
	}
	lines(beta.average.quantiles[1:2,v],rep(k+2,2),col="grey",lwd=2)
	points(beta.average.est[v],k+2,pch=16,cex=1)
	
	# Chief Justice Divisions #
	text(3,1.5,"Vinson",pos=2,cex=0.75)
	abline(h=2.5)
	text(3,8,"Warren",pos=2,cex=0.75)
	abline(h=13.5)
	text(3,17,"Burger",pos=2,cex=0.75)
	abline(h=20.5)
	text(3,24,"Rehnquist",pos=2,cex=0.75)
	abline(h=27.5)
	text(3,29,"Average",pos=2,cex=0.75)
	box("figure",col="grey")	
	

}
	box("outer",col="grey")	
dev.off()






















###############################
## Run JAGS Models for Joins ##
###############################



keepCases <- ifelse(rowSums(is.na(AuthorDistance)) + rowSums(is.na(MedianDistance)) + rowSums(is.na(MajMedianDistance)) + rowSums(is.na(ChiefDistance))==0,1,0) & 
	ChiefInMajority & 
	(rowSums(JoinMatrix, na.rm=TRUE)>0) & # someone joins
	(apply(latent.preferences,1,sd) < 2) # removes disjoint cases with preference estimates that are not identified relative to the others



# load data and run simulations
forJags <- list(
  Join = JoinMatrix[keepCases,],
  Author = Author[keepCases],
  MedianDistanceEst = MedianDistance[keepCases,],
  MajMedianDistanceEst = MajMedianDistance[keepCases,],
  CutDistanceEst = CutDistance[keepCases,],	
  ChiefDistanceEst = ChiefDistance[keepCases,],
  AuthorDistanceEst = AuthorDistance[keepCases,],
  MedianDistancePrec = MedianDistanceSE[keepCases,]^(-2),
  MajMedianDistancePrec = MajMedianDistanceSE[keepCases,]^(-2),
  CutDistancePrec = CutDistanceSE[keepCases,]^(-2),	
  ChiefDistancePrec = ChiefDistanceSE[keepCases,]^(-2),
  AuthorDistancePrec = AuthorDistanceSE[keepCases,]^(-2),
  NotInMajority = NotInMajority[keepCases,],
  MajorityVotes = MajorityVotes[keepCases],
  NaturalCourt = as.numeric(voteDetails$NaturalCourt)[keepCases],
  ChiefOfNaturalCourt = ChiefOfNaturalCourt,
  Ncases = table(keepCases)[2],
  Njustices = dim(MedianDistance)[2],
  Ncourts = length(levels(voteDetails$NaturalCourt)),
  Nchiefs = length(unique(ChiefOfNaturalCourt))
)



ModelNames <- c("Main","1","2","3","4","5")  

if (RunMCMC){
  for (current.model in ModelNames){
    
    current.model.string <- paste("OpinionJoinMeasurementError", current.model,".jags",sep="")
    
    model.jags <- jags.model(file=current.model.string
                             ,data=forJags
                             ,n.chains=MCMCChains,n.adapt=round(MCMCBurn/4))
    
    update(model.jags, round(3*MCMCBurn/4))
    
    posterior <- jags.samples(model.jags, 
                              c("beta","beta2","beta3"
                                ,"alpha","alpha2","alpha3"
                                ,"sigmaalpha","sigmaalpha2"
                                ,"sigmabeta","sigmabeta2"
                                ,"sigmagamma","sigmagamma2"
                                ,"deviance","pD"
                              ), 
                              n.iter=MCMCSample, thin = MCMCThin, progress.bar="text")
    
    current.posterior.string <- paste("Output/OpinionJoinMeasurementError", current.model,".RData",sep="")
    
    save(posterior, file=current.posterior.string)	
    
    DIC <- mean(posterior$deviance) + mean(posterior$pD)	
    print(paste("Model",current.model,"DIC:",DIC))
    
  }	
}


########################
## Analyze Full Model ##
########################	

current.posterior.string <- paste("Output/OpinionJoinMeasurementErrorMain.RData",sep="")
load(file=current.posterior.string)		

DIC <- mean(posterior$deviance) + mean(posterior$pD)	
print(paste("DIC:",DIC))

alpha3.est <- apply(posterior$alpha3,1,mean)

alpha2.est <- apply(posterior$alpha2,c(1,2),mean)
rownames(alpha2.est) <- colnames(Y)
colnames(alpha2.est) <- colnames(Y)

alpha.est <- apply(posterior$alpha,c(1,2,3),mean)
dimnames(alpha.est)[[1]] <- dimnames(alpha.est)[[2]] <- colnames(Y)
dimnames(alpha.est)[[3]] <- as.character(unique(voteDetails$NaturalCourt))

beta3.est <- apply(posterior$beta3,1,mean)
beta3.quantiles <- t(apply(posterior$beta3,1,quantile,c(0.025,0.975)))
print(cbind(beta3.est,beta3.quantiles))

beta2.est <- apply(posterior$beta2,c(1,2),mean)
rownames(beta2.est) <- c("Distance to Cutpoint","Distance to Court Median","Distance to Majority Median","Distance to Chief Justice","Distance to Author")
colnames(beta2.est) <- as.character(unique(voteDetails$Chief))
#beta2.quantiles <- t(apply(posterior$beta2,c(1,2),quantile,c(0.025,0.975)))
print(beta2.est )

beta.est <- apply(posterior$beta,c(1,2),mean)
beta.quantiles <- apply(posterior$beta,c(1,2),quantile,c(0.025,0.975))
rownames(beta.est) <- c("Distance to Cutpoint","Distance to Court Median","Distance to Majority Median","Distance to Chief Justice","Distance to Author")
colnames(beta.est) <- as.character(unique(voteDetails$NaturalCourt))
print(beta.est)



#####################
## Plot Full Model ##
#####################	

pdf("Output/JoinCoefficientPlotMultivariate.pdf",12.5,4)
label.side <- c(-1,-1,1,1,1)
n.periods <- dim(beta.est)[2]
par(mfrow=c(1,5),mar=c(4,4,4,4))
for (v in 1:5){
	plot(0,0,xlim=c(-5,5),ylim=c(0,n.periods+4),type="n",axes=FALSE,xlab="Coefficient",main=rownames(beta.est)[v],ylab="Natural Court")
	axis(1)
	abline(v=0)
	for (k in 1:n.periods){
		lines(beta.quantiles[1:2,v,k],rep(k,2),col="grey")
		points(beta.est[v,k],k,pch=16,cex=0.75)
	}
	lines(beta3.quantiles[v,1:2],rep(k+2,2),col="grey",lwd=2)
	points(beta3.est[v],k+2,pch=16,cex=1)
	
	# Chief Justice Divisions #
	text(label.side[v]*5,1.5,"Vinson",pos=3-label.side[v],cex=0.75)
	abline(h=2.5)
	text(label.side[v]*5,8,"Warren",pos=3-label.side[v],cex=0.75)
	abline(h=13.5)
	text(label.side[v]*5,17,"Burger",pos=3-label.side[v],cex=0.75)
	abline(h=20.5)
	text(label.side[v]*5,24,"Rehnquist",pos=3-label.side[v],cex=0.75)
	abline(h=27.5)
	text(label.side[v]*5,29,"Average",pos=3-label.side[v],cex=0.75)
	
}
box("outer",col="grey")
dev.off()







pdf("Output/JoinCoefficientPlotUnivariate.pdf",12.5,4)
label.side <- c(-1,1,1,1,1)
n.periods <- dim(beta.est)[2]
par(mfrow=c(1,5),mar=c(4,4,4,4))
for (v in 1:5){
	current.posterior.string <- paste("Output/OpinionJoinMeasurementError", v,".RData",sep="")
	load(file=current.posterior.string)	
	
	DIC <- mean(posterior$deviance) + mean(posterior$pD)	
	print(paste("DIC:",DIC))
	
	alpha3.est <- apply(posterior$alpha3,1,mean)
	
	alpha2.est <- apply(posterior$alpha2,c(1,2),mean)
	rownames(alpha2.est) <- colnames(Y)
	colnames(alpha2.est) <- colnames(Y)
	
	alpha.est <- apply(posterior$alpha,c(1,2,3),mean)
	dimnames(alpha.est)[[1]] <- dimnames(alpha.est)[[2]] <- colnames(Y)
	dimnames(alpha.est)[[3]] <- as.character(unique(voteDetails$NaturalCourt))
	
	beta3.est <- apply(posterior$beta3,1,mean)
	beta3.quantiles <- t(apply(posterior$beta3,1,quantile,c(0.025,0.975)))
	print(cbind(beta3.est,beta3.quantiles))
	
	beta2.est <- apply(posterior$beta2,c(1,2),mean)
	rownames(beta2.est) <- c("Distance to Cutpoint","Distance to Court Median","Distance to Majority Median","Distance to Chief Justice","Distance to Author")
	colnames(beta2.est) <- as.character(unique(voteDetails$Chief))
	#beta2.quantiles <- t(apply(posterior$beta2,c(1,2),quantile,c(0.025,0.975)))
	print(beta2.est )
	
	beta.est <- apply(posterior$beta,c(1,2),mean)
	beta.quantiles <- apply(posterior$beta,c(1,2),quantile,c(0.025,0.975))
	rownames(beta.est) <- c("Distance to Cutpoint","Distance to Court Median","Distance to Majority Median","Distance to Chief Justice","Distance to Author")
	colnames(beta.est) <- as.character(unique(voteDetails$NaturalCourt))
	print(beta.est)

	# Calculateaverage of natural-court specific coefficients
	beta.average.chain <- apply(posterior$beta,c(1,3),mean)
	beta.average.est <- apply(beta.average.chain,1,mean)
	beta.average.quantiles <- apply(beta.average.chain,1,quantile,c(0.025,0.975))

	
	plot(0,0,xlim=c(-5,5),ylim=c(0,n.periods+4),type="n",axes=FALSE,xlab="Coefficient",main=rownames(beta.est)[v],ylab="Natural Court")
	axis(1)
	abline(v=0)
	for (k in 1:n.periods){
		lines(beta.quantiles[1:2,v,k],rep(k,2),col="grey",lwd=1)
		points(beta.est[v,k],k,pch=16,cex=0.75)
	}
	lines(beta.average.quantiles[1:2,v],rep(k+2,2),col="grey",lwd=2)
	points(beta.average.est[v],k+2,pch=16,cex=1)
	
	# Chief Justice Divisions #
	text(label.side[v]*5,1.5,"Vinson",pos=3-label.side[v],cex=0.75)
	abline(h=2.5)
	text(label.side[v]*5,8,"Warren",pos=3-label.side[v],cex=0.75)
	abline(h=13.5)
	text(label.side[v]*5,17,"Burger",pos=3-label.side[v],cex=0.75)
	abline(h=20.5)
	text(label.side[v]*5,24,"Rehnquist",pos=3-label.side[v],cex=0.75)
	abline(h=27.5)
	text(label.side[v]*5,29,"Average",pos=3-label.side[v],cex=0.75)
	box("figure",col="grey")	
}
dev.off()























